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1.1 Introduction 

This chapter is devoted to the computation of cquihbrium (thermodynamic) 
properties of quantum systems. In particular, we will be interested in the sit- 
uation where the interaction between particles is so strong that it cannot be 
treated as a small perturbation. For weakly coupled systems many efficient 
theoretical and computational techniques do exist. However, for strongly inter- 
acting systems such as nonideal gases or plasmas, strongly correlated electrons 
and so on, perturbation methods fail and alternative approaches are needed. 
Among them, an extremely successful one is the Monte Carlo (MC) method 
which we are going to consider in this chapter. 

1.2 Idea of Path integral Monte Carlo 

If we perform classical simulations of a system in equilibrium, we usually 
need a Boltzmann-type probability distribution, pB ~ e~^^"''-^'> /Z , (/? = 
l/ksT) and then the Monte Carlo method [1, 7] can be used to sample the 
particle coordinates R= (ri, r2, . . . , rjv). Now the question arises what is the 
appropriate probability density in the quantum case. The answer is provided 
by the density operator p. Consider a general expression for thermodynamic 
averages in statistical thermodynamics. The A'^— particle density matrix p{(3) 
contains the complete information about the system with the observables 
given by {Tr denotes the trace) 




(1.1) 
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This expression is simplified if the physical observable O is diagonal in the 
coordinate representation, i.e. {B!\6 p{P)\R) ^ {R\d p{l3)\R) 5{R' - R). 

As in the classical case we need to perform an integration over 3A^ (or 
more) degrees of freedom. The problem here, compared to the classical case, 
is that, in general, we do know an analytical expression for the A^-particlc 
density operator in the coordinate representation which can be substituted in 
Eq. (1.1). This problem was first overcome by Feynman [2]. The key idea is 
to express the unknown density operator, 

{R\p{l3)\R') with p = e-^", (1.2) 

at a given inverse temperature /? by its high-temperature asymptotic which is 
known analytically. But the price to pay is high: instead of an already com- 
plicated 3 A'^— dimensional integral, now it expands to much higher dimensions 
(3iV- M), where M is an integer which in practice is chosen as 1 < M < 3000. 

A. Group property of DM 

One simple and straightforward strategy is to use the group property of the 
density matrix. It allows to express the density matrix at low temperatures 
in terms of its values at higher temperature, i.e. 

p{R,R';Pi+P2) = {R\e-'-'^'+'^'^"\R') = jdRi {Rle-"'" \Ri) {Ri\e-^'^ \R') 

= JdRip{R,Ri;pi)p{RuR']/32). (1.3) 

The next step is a generalization: use the group property M times 

p = exp[-pH] = exp[-ApH] . . . exp[-Aj3H], Ap = p/M. (1.4) 

This means that the density operator p is expressed as a product of M new 
density operators, e~^^^ , each corresponding to an M times higher temper- 
ature. 

Finally, using Eq. (1.4) we can write^ the off-diagonal matrix element as 
(for any fixed end-points R and R') 

piR,R';p) = jdRidR2...dRM-i 

p{R, Ri; A(3)p{RiR2; Ap) . . . p{Rm-i,R'; A(3), (1.5) 

where M factors are connected by M — 1 intermediate integrations. 

B. High temperature approximation 

The Eqs. (1.4) and (1.5) are exact at any finite M as long as we use exact ex- 
pressions for the high-temperature A^-particle density matrices, p{Ri-\,Ri; Ap) 

^ The total dimension of the integral, (M — 1) SAT, may be very large. The success 
of the method relies on highly efficient Monte Carlo integration. 
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Unfortunatly, they are unknown, and to proceed further we need to introduce 

approximations . 

Our approximation will be based on Trotter's theorem (1959) applied to a 
general Hamiltonian, H = T + V, which contains both kinetic and potential 
energy operators, i.e. 



p = lim 

M^oo 



^-AfST^-A/3V 



M 



^-Al3T^-A(iV 



M 



g-A/3Tg-Zi/3V 



M 



(1.6) 

+ 0(1/M). 
(1.7) 



Note that T and V do not commute giving rise to the commutator, [T, T^], 
which is only the first term of a series^. Neglecting the terms [T, V^] gives 
an error of the order O [1/M]. This error can be made arbitrarily small by 
choosing a sufficiently large number of factors M. 

Using the Trotter result (1.7), we immediately obtain an approximation 
for high temperatures^ 



= A^3^exp 



-^{Ri- Ri+if - APV{Ri; Ap) 



(1.8) 



where = {2iTh^A(3/m)2 is the De Broglie wavelength. Substituting 
Eq. (1.8) in Eq. (1.5) we get our final result for low-temperatures 



p{R, R'; p) = jdRi... dRM-i e 



-{R,-Ri+i) 



E AI3V{R,) 

i=0 



(1.9) 



with the boundary conditions: Rq = R and Rm = R' ■ Hence, we have suc- 
ceeded in the analytical representation of the A^-particle density matrix. This 
representation is then used in the numerical evaluation of this integral with 
the help of the Monte Carlo based algorithm. 

C. Visualization of diagonal elements of the density matrix 
As we can see from Eqs. (1.5) and (1.9), now all A'^ particles have their own 
images on M different planes (or "time slices"). We can view these images 
(for each particle 3 ■ M sets of coordinates) as a "trajectory" or a "path" 
in the configurational space. The inverse temperature argument (3 can be 
considered as an imaginary time of the path. The set of M time slices is 
ordered along the /3-axis and separated by intervals A[3. On Fig. 1.1 we show 
typical configurations of particle trajectories which contribute to the diagonal 
density matrix element, Eq. (1.5) with R = R' . The full DM p{R,R;f3) is 



^ Double, triple and high-order commutators in a series have as a prefactor higher 

powers A/3" and they can be dropped in the limit A/S —^ 0. 
^ Other more accurate high temperature approximations are discussed in Ref. [6, 7]. 
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Fig. 1.1. Live picture from a computer simulation. Five particles in a 2D parabolic 
potential. Each particle is presented as a continuous path (obtained by a smooth 
interpolation through a discrete set of M = 100 points). Right fig. shows how the 
paths are stretched along the time (/3) axis. Particles 1 and 2 are in a pair exchange. 
Labels show particle indices. 



obtained after integration over all possible path configurations with the fixed 
end points {R = R'). 

If we look at the final analytical result for the high temperature density 
matrix, Eq. (1.9), we recognize the usual Boltzmann factor with some effective 
action in the exponent. This action describes two types of interaction. The 
first term, 

N M-1 

^Y^H^i-^i+if (1-10) 

comes from the kinetic energy density matrices of free particles (j denotes 
the particle index). This energy can be interpreted as the energy of a spring. 
Us = |(Z\a;)^. Changing one of the coordinates is equivalent to a change of 
the spring energy. These springs provide that the nearest points on the path 
are usually at some average distance proportional to \a- With temperature 
reduction the average size of the path increases with \a- 

The second term Al3V{Ri) in Eq. (1.9) adds interactions to the system 
(e.g. an external potential or inter-particle pair interaction) 

M~l I N M-1 M-1 \ 

J2 ^PV{R,) - E E ^-*(^') + E E ^p-(r^rf) . (1.11) 

Each potential term depends only on the particle coordinates on the same 
time slice, i.e. (r^, , . . . , y^). As a result the number of pair interactions on 
each time slice, N[N — l)/2, is conserved. 



E 



R, 



-1? = ■ 
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In all expressions above we have considered particles as distinguishable 

ones. Generalization to quantum particles obeying Bosc statistics is considered 
in Sec. 1.4. The case of Fermi statistics is discussed e.g. in Ref. [4, 5, 6, 7]. 

1.3 Basic numerical issues of PIMC 

Having the general idea of the PIMC simulations we are ready to formulate 
the first list of important issues which we need to solve: 

A. How to sample the paths 
It is necessary to explore the whole coordinate space for each intermediate 
point. This is very time consuming. To speed up convergence: move several 
slices (points of path) at once. 

Key point: sample a path using mid-points Rm and a consequent iteration 
{bisection), see Fig. 1.2(b). 



Using the definition: < i < /3, r = ioAp [io = 1,2,3, ...], R = R{t), R' = 
R{t + 2r), Rm = R{t + r); Guiding rule to sample a mid-point Rm is 



In practice, we can neglect in the sampling distribution the potential energy 
and use only the ratio of the free-particle density matrices. As a result we 
get a Gaussian distribution with the mean, R = {R + R')/2 and the variance, 
fT^ = f{^T/2m. This will lead to 100% acceptance of sampling for ideal systems 
(and close to one for a weakly interacting system). 

For strongly interacting systems the overlap of the paths sampled from 
the free-particle distribution (1.12) results in large increase of the interaction 
energy and in a poor acceptance probability at the last level of the bisection 
sampling [6, 7]. This can be improved by using the optimized mean and the 
variance 



R={R + R')/2 + ar^^^, al=h^T/2m+(—) AV{R), (1.13) 



which account for local interaction strength (nearest neighbor interaction). 
Advantages of the bisection sampling method [6, 7]: 

• Detailed balance property is satisfied at each level. 

• We do not waste time on moves for which paths come close and the po- 
tential energy strongly increases (for repulsive interaction). Such configu- 
rations are rejected already on early steps. 

• Computer time is spent more efficiently because we consider mainly con- 
figurations where the acceptance rate is high. 

• Sampling of particle permutations is easy to perform. 



P{Rm) = 



{R\e-^"\Rm){Rm\e- 
{R\e-^^"\R') 



\R') 



(27ra2)-''/2g-(K„-ii)V2.J_ (1_12) 
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B. Choose action as accurate as possible 

It is better to use effective interaction potentials which take into account 
two, three and higher order correlation effects. More accurate actions help to 
reduce the number of time slices by a factor of 10 or more. Prom Eq. (1.9) we 
can note, that the first free-particle terms can be considered as a weight over 
all possible random walks (Brownian random walks) in the imaginary time (3 
with the ends points R and R'. In the limit M — > oo we directly obtain the 
Feynman-Kac relation 

/ - J dtV{R.{t))\ 

p{R,R';P) = po{R,R';P)(e « ) . (1.14) 

\ / FK 

In the quasi-classical limit (/3 — > 0), only the classical path is important, 
Ro{t) = (l — R+ ^R' , which leads to the semi-classical approximation of 
the high-temperature density matrix 

- f dtV{Ro{t)) 

p{R,R';Ap)^poiR,R';Ap)e o , (1.15) 

which is already much better compared to Eq. (1.8) with the substitution of 
classical (in many cases divergent) potentials. 

For systems with pair interactions, in the limit of small Z\/3, the full density 
matrix (1.14) can be approximated by a product of pair density matrices 

/ -JdtV{R{t))\ I -/dtyp„„[r,(t),r,(t)]\ 

\ / FK j<k \ I FK 

This is used in practice as the pair approximation. It supposes that on the 
small time interval the correlations of two particles become independent 
from the surroundings. Different derivations of the effective pair potential 
(average on the l.h.s of Eq. (1.16)) have been proposed in literature [8, 9, 10]. 

Implementation of the periodic boundary conditions leads to further mod- 
ifications, see e.g. Refs. [11, 12, 13]. 

C. How to calculate physical properties 

There are different approaches for calculating expectation values of physical 
observables, such as the energy, momentum distribution, etc., which are called 
estimators. In each particular case convergence can be improved by choosing 
the proper estimator. Consider, for example, the thermodynamical estimator 
of the internal energy 

a ^. IdZ 1 f-dp{R,R;P) ,^ ^ , 

After direct substitution of Eq. (1.9) for p, one obtains^ 

* Tliis is only valid for particles with Boltzmann statistics. For fermions one has 
to include additional terms related to the /3-derivative of the exchange determi- 
nant [16, 7]. 
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»=0 / p(fi:/3) \ i=0 I p(R;l3) 

(1.18) 

This form of the energy estimator has a much larger statistical variance ag 
compared to the virial estimator [14]. Since the statistical error in Monte Carlo 
simulations decreases as SE « (Js/V^mc (with Nmc being the number of 
MC-steps), with the direct estimator (1.18) one usually needs several times 
more MC runs to get the same accuracy as given by the virial estimator. 

One of the approaches to obtain the virial estimator of energy is to intro- 
duce temperature dependent coordinates [15], i.e. Ri = Rq+Xa X^to=i * ~ 
1, . . . , M — 1. Here is a set of unit vectors, and Ro is a set of particle coordi- 
nates at the zero time slice {Rq = R). Once this has been done, the estimator 
takes the form [7] 

One can note at once, that for weakly interacting systems at high tempera- 
tures, the virial result (1.19) directly gives the classical kinetic energy (first 
term) and does not depend on the chosen number of time slices M, whereas 
in the direct estimator (1.18) we get this result by calculating the difference 
of two large terms which are diverging as M — > oo. 

D. Acceptance Ratio 

When we try different kinds of moves in the Metropolis algorithm, it can turn 
out that some moves will be frequently rejected or accepted. In both cases, we 
loose the efficiency of the algorithm. For a sufficiently long time (number of 
MC steps) we will be trapped in some local region of phase space and will not 
explore the whole space in a reasonable computer time. In practice, the pa- 
rameters of the moves are usually chosen to get an acceptance ratio of roughly 
50%. For different kinds of moves in PIMC (particle displacement, path de- 
formation, permutation sampling) such an acceptance ratio is preferable and 
construction of good apriori sampling distributions is crucial. 

A detailed discussion of these topics is beyond the scope of the present 
paper and is covered in other publications, e.g. Rcfs. [6, 7]. 

E. Quantum exchange. PIMC for bosons/fermions 

Now we come to "real" quantum particles. As we have already discussed, 
properties of a system of A'^ particles at a finite temperature T arc determined 
by the density operator. Due to the Fermi/Bose statistics the total density 
matrix should be (anti) symmetric under arbitrary exchanges of identical par- 
ticles (e.g. electrons, holes, with the same spin projection), i.e. we have to 
replace p p^^^ for fermions/bosons. As a result the full density matrix will 
be a superposition of all A^! permutations of A'' identical particles. Consider 
the case of two types (e,h) of particles with numbers A^e, 



8 Alexei Filinov, Jens Boning, and Michael Bonitz 

(1.20) 

where Pe(h) is the parity of a permutation (number of equivalent pair transpo- 
sitions) and f'e(h) the permutation operator. We directly see that, for bosons 
all terms have a positive sign, while for fcrmions the sign of the prefactor 
alternates depending whether the permutation is even or odd. 

In the last case a severe problem arises. The Metropolis algorithm gives 
the same distribution of permutations for both Fermi and Bosc systems. The 
reason is that, for sampling permutations, we use the modulus of the off- 
diagonal density matrix, |p(ii. Pi?; /3)|. 

• Bose systems: all permutations contribute with the same (positive) sign. 

Hence with the increase of the permutation statistics, accuracy in the 
calculation of the density matrix increases proportionally. 

• Fermions: essential cancellation of positive and negative terms (corre- 
sponding to even and odd permutations), both are close in their absolute 
value. Accurate calculation of this small difference is drastically hampered 
with the increase of quantum degeneracy (low T, high density). The con- 
sequences are large fluctuations in the computed averages. This is known 
as the fermion sign problem. It was shown [5] that the efhciency of the 
straightforward calculations scales like e"^^'^^^, where AF is the free en- 
ergy difference per particle of the same fermionic and bosonic system, and 
N is the particle number. 

F. Numerical sampling of permutations 

Fermi and Bose statistics require sampling of permutations, sec Eq. (1.20), in 
addition to the integrations in real space. From the A^! possibilities, we need 
to pick up a permutation which has a non-zero probability for a given particle 
configuration. 

To realize a permutation we pick up, along the /3-axis, two end-points 
{Ri, Ri+io} with io = 2'^^ {I = 1,2,. . .). Although the permutation operator 
P in Eq. (1.20) acts on the last time-slice, Re(h)^ the permutation of the paths, 
{Ri, PRi^i„} can be carried out at any time slice because 
the operator P commutes with the Hamiltonian. In a permutation move {k 
permuted particles) the path coordinates between the fixed points Ri and 
Ri+i„ are removed and new paths connecting one particle to another (new 
k links) or a new path connecting a particle on itself (if a given particle 
undergoes the identity permutation) are sampled. 

It is evident, that a local permutation move consisting of a cyclic exchange 
of > 2 neighboring particles will be more probable than a global exchange 
involving a macroscopic number of particles, and, in general, the probability 
of exchange will decrease with the increase of k. The most probable are local 
updates: permutations of only few (2,3,4) particles. Moreover, any of the N\ 
permutations can be decomposed in a sequence of successive pair transpo- 
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sitions (two particle exchange), and we can explore the whole permutation 

space by making only local updates which have a high acceptance ratio. 
In MC simulations we choose as the sampling probability of permutations 

T{p p') = PkinjRi, PRi+i„;iQ^P) 
where we have used the product of the k one-particle density matrices 

Pk^n{R^,PR^+^o,■,^oAP) (X CXp I " ^ 7r(r^' - Pvi+J^ / (ioXl) • (1-22) 

Here denotes the neighborhood of the current permutation P from which 

the permutation P' is sampled. For example, for the two-particle exchange, 
J7(P) equals to number of neighbors of the given particle in the range of 
several De Broglie wavelengths, X{t), t = io^P, which are possible candidates 
for the exchange. 

To satisfy the detailed balance principle, we make a final decision^ about 
the sampled permutation using the acceptance probability 



A{P P') = min 



/ ^ Pkin 

{Ri, PRi+ig-jioAfi) 

E Pkin {Ri , P'Ri+io ; ioAp) 

p'en(P') 



(1.23) 



where fi{P') is the neighborhood of the new permutation P' . If the neighbor- 
hoods of the current and new permutation are equal, the acceptance proba- 
bility is one. 

As an illustration, in Fig. 1.2 we show a world line picture of five particles. 
Particle indices in Fig. 1.2 (a) and (c) are placed near the starting and end 
point of the particle trajectories. Hence, when the sequence of indices at m = 
and m = 100 does not coincide, sec Figs. 1.2 (a)(b), the particles are permuted. 

As we can see from Fig. 1.2(a) the paths of particles "1" and "2" are closed 
(two identical permutations), three other particles are in one cyclic exchange, 
and the whole permutation can be denoted as {1, 2. 5, 3. 4} (as we can see the 
end of the path "3" coincides with the beginning of path "5" , the end of path 
"4" coincides with the beginning of path "3" and the path "5" ends up at the 
starting position of path "4" . Now we decide to make a transposition between 
particles "1" and "4" . To do this we choose randomly time slices where new 
paths will be sampled. In our case it was m = 17—33. First, we exchange the 
edge points at the time slice m = 33, i.e. rf = and r4^'^ = rf' . Hence the 
position of the edge points is not sampled, they are a part of the unchanged 
trajectories. Once the initial and final points are chosen, we use the bisection 
algorithm to sample two paths connecting edge points, see Fig. 1.2(b). 



® The sampled permutation can be rejected earlier when the new paths connecting 
Ri and PRi+i^ are sampled with the bisection algorithm [6, 7]. 
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Fig. 1.2. (a),(c) The Y-coordinates of five identical particles as a function of the 
time-slice number m. Labels show particle indices. Thick gray and light gray lines 
show the paths of the particles "1" and "4" which are exchanged by sampling new 
paths at time-slices m = 17 — 33 (these time-slices are in the region between two 
dashed lines), (b) Sampling of new paths using the bisection algorithm for the elec- 
trons "1" and "4" . The new paths are constrained at the time-slices m = VJ — 33. 
Old (new) paths are shown by lines (circles). The filled circles show two mid-points 
sampled at the level I — 1 (center of the interval, m — 25) and four other mid-points 
for sub-intervals [17, 25] and [25, 33] sampled at level I — 2. Open circles show final 
new paths for two particles obtained with the sampling at levels Z = 3, 4 and the 
transposition, i.e. by exchanging the paths starting from m = 33 up to the end 
point, m = 100. 

1.4 PIIVIC for degenerate Bose systems 

Currently much experimental activity is devoted to the study of ensembles 
of dilute gases of Bose atoms and optically excited indirect excitons in sin- 
gle/double well nanostructures, e.g. Refs. [17, 18] and references therein. The 
most exciting is certainly the possibility to observe signatures of Bose conden- 
sation and superfluidity. The essential point of these experiments is that the 
number of trapped atoms are limited to a few ten thousand particles and one 
should expect significant deviation (finite size effects) from the macroscopic 
limit, leading e.g. to a "softening" of the condensate fraction curve in the 
transition region and also to a shift of the critical temperature to lower val- 
ues. This is particularly important for the case of a few hundreds of particles 
which become accessible to a direct theoretical investigation using Quantum 
Monte Carlo approaches which allows to treat many-body correlation effects 
from "first principles". 
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In PIMC, as was shown by Feynman [2], the Bose statistics manifest itself 
as a special topology of the particle trajectories which can form macroscop- 
ically large permutation cycles. The free external parameters, like temper- 
ature, density, interaction strength, have a direct influence on these cycle 
distributions and, hence, on the superfluid and condensate fractions. Below 
we demonstrate how the latter can be easily related to the statistics of path 
configurations sampled by PIMC. 

To be more specific in the discussion below we consider a system of trapped 
bosons with Coulomb interaction described by the Hamiltonian (which can 
be also reduced to the dimensionless form, H, using the system of units: 
r ^ f = {r/aB),E E = {E/Ra) with as = h^e/me^ and Ha = e^/eaB) 

* = A. + E7^. «. = i: -|sVj. + f.^.?), (1.24, 
iJ^H/Ha^f:l(-Vi, + ln»)+f:±. (1.25, 

i=l ^ ' i<j 

For the mesoscopic trapped system as considered here the density is controlled 
by the harmonic trap frequency lu and is characterized by the coupling pa- 
rameter A = {e^ / elo) / {hw) = lo/as with Iq = h/mco. In the limit (A —^ oo) 
the external potential vanishes, while for (A — > 0) the Coulomb interaction 
can be neglected (formal transition to non- interacting bosons). 
A. Superfluidity 

First, we consider the fraction of the superfluid (mass) density 7^ = ps/p 
which, within the Landau two-fluid model is computed from the classical and 
quantum momenta of inertia, Ic and Ig, according to 7s = 1 — Iq/Ic [19]. The 
quantities Ic and Iq are eS'ectively computed in PIMC simulations [20] from 
the area enclosed by the particle paths A, using 

(1.26) 

where A'' is the particle number, M the number of time slices used in the 
path integral presentation and (. . .) denotes the thermal averaging with the 
symmetric A'^— particle diagonal density matrix, i.e. 

(•••) = |y J dv,dv2...dvM (...) p^{vur2,...,VM;P). (1.27) 

This formula has been derived [20] for finite systems by assuming that particles 
are placed in an external field, e.g. a rotating cylinder. Then one considers 
that the system is put in a permanent slow rotation with the result that 
the normal component follows the imposed rotation while the superfluid part 
stays at rest. The effective moment of inertia is defined as the work done to 
rotate the system by a unit angle. 



12 



Alexei Filinov, Jens Boning, and Michael Bonitz 



For macroscopic systems the path area formula (1.26) can be modified [6, 
21]. Instead of a filled cylinder, one considers two cylinders with the radius 
R and spacing d, where d <^ R. Such a torus is topologically equivalent to 
the usual periodic boundary conditions. As a result we have: Ic = mNR"^ and 
Az = WR/2, where W is the winding number, defined as the flux of paths 
winding around the torus and crossing any plane 

p -2A/3iV' ^-^io 

B. Off-diagonal long-range order 
The magnitude of off-diagonal long-range order is, in macroscopic systems, 
also directly accessible with PIMC. It is characterized by the asymptotic be- 
havior of the single-particle off-diagonal density matrix 

p{vi,r'^;(3) = V/Z j dra . . . drjv /(ri, ra, . . . , ta,, r'l, rs, . . . , tat; /3), (1.29) 
no(/3) = lim p(ri,r;;/3), (1.30) 

where no is the fraction of particles in the condensate and V is the volume of 

the simulation cell. For a homogeneous isotropic system, p(ri, v'l) = p{\ri—r[\) 
and, by taking the Fourier transform of an off-diagonal element, one obtains 
the momentum distribution 

P^^'^^jirr I <ri-r;)e-''(-^-^i)p(|ri-r;|;/3), (1.31) 

which shows a sharp increase at zero momentum when the temperature drops 
below the critical temperature Tc of Bose condensation. 

Obviously, a finite trapped system of particles considered in real experi- 
ments behaves differently. The radial density is strongly inhomogeneous with 
the highest value at the trap center. However, these systems do represent 
an analog of the homogeneous macroscopic system in the angular direction 
(for traps with angular symmetry as in the case (1.24)). Hence, the macro- 
scopic formulas (1.30) and (1.31) should be modified in an appropriate way 
and the corresponding momentum distribution, the condensate fraction and 
superfluidity get an additional dependence on the radial distance from the 
trap center. 

As follows from Eq. (1.29), for the numerical evaluation of the single- 
particle density matrix one should allow that one of the N simulated particles 
has an open path, e.g. ri ^ r^. The paths of the other N — 1 particles can 
close on their beginning (identical permutation) or at the start of another 
particle's path. The coordinates ri , r[ are independent but their probability 
is given by the TV-particle density matrix. In simulations we record a histogram 
(distribution) given by 



drjjt) 
dt 



(1.28) 
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p(r,r';/?)cx(5(ri-r)5(r;-r'))^, 
W = p^(ri, r2, . . . , rjv, r'l, r2, . . . , tn; P)/Z', 



(1.32) 
(1.33) 



(Z' is the normalization factor) which is then used to obtain the momentum 
distribution (1.31). The probabiUty W is sampled using the path integral 
representation of p^' . 

Recently a new method to sample the single-particle density matrix (1.29), 
by generalization of the conventional PIMC to the grand canonical ensemble, 
has been proposed [22]. A worm algorithm allows for a simultaneous sam- 
pling of both diagonal configurations contributing to the partition function 
and off-diagonal ones which contribute to the one-particle Matsubara Green 
function. The method has been recently applied to study of Bose condensa- 
tion in crystalline '^He and superfluidity in para- hydrogen droplets [23] , where 
high efficiency in sampling of long permutation cycles (practically unaffected 
by system size) and significantly improved convergence in the calculation of 
superfiuid properties has been demonstrated. 

1.5 Numerical results 

A. Reference case: Ideal system 

For the noninteracting Bose gas the total Hamiltonian is a sum over inde- 
pendent single particle Hamiltonians, H = Hg, Eq. (1.24). If we then calcu- 
late the partition function as the trace over the symmetrized density matrix 
p^{R, R; P), it can be written as a product of N independent integrals 



where each integral contains only the single-particle density matrix. 

Any permutation can be broken into exchange cycles. A given permutation 
consists of a set of cycles {Ci, C2, . . . , Cat}, where Cq denotes the number of 
cycles with length q. As the particles are unlabeled, different permutations 
resulting in identical cycle configurations give identical contributions to the 
partition function. We can thus rewrite the sum over permutations as a sum 
over cycle configurations with an additional factor taking into account the 
number of possible realizations. By applying the group property of the density 
matrix the integrals over the single-particle density matrices within a cycle of 
length q are reduced to a single integral over the density matrix at the g-times 
lower temperature with the result [3] 




(1.34) 




(1.35) 



restr. 
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with the restriction J2q^i qCq = N, for the possible cycle configurations. 

The same procedure can be applied for the calculation of any thermal 
average. In particular, for the mean cycle occupation number (C,) we obtain 

restr. 

Plugging this formula into the constraint on the possible cycle configurations 
leads to a powerful recursion relation for the partition function [4, 24] 



^^('^) = ^ E ZN-MZim- (1-37) 
For the energy of the system we obtain (r = 1//3) 

{E) = r^^MEM. ilCq) . (1.38) 

9=1 

The occupation number Ni of an arbitrary energy level Ei can be cal- 
culated as the derivative of the partition function with respect to l3Ei, see 
Ref. [24] 

1 dZ^m _ ^ e-^^^' 

z^)-dm)-^,z^)^'^'^- ^'-''^ 

This formula yields the number of particles in the condensate, when taken at 
the ground state energy Eq. 

Within linear response theory, the quantum mechanical moment of inertia 
7q is defined as the system response to rotations. A careful analysis for non- 
rotating situations with (Lz) = of a system confined in a 2d spherical 
harmonic trap yields [25] 



AT 



I, = (3 {Ll) = 2h- Y: n'-l-q,n.^2 ilCq) . (1.40) 

9=1 ^ ^ 

The classical moment of inertia 7c of the same system becomes 

j=l 9=1 

The left side of Fig. 1.3 shows a typical cycle configuration, i.e. the prob- 
ability, that a particle belongs to an exchange cycle of length q. Clearly, ana- 
lytical and numerical results are in good agreement. Both indicate that every 
possible permutation is equally probable at low temperatures, whereas only 



q0hLJ 
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Temperature, kT/hra Temperature, kT/hm 



Fig. 1.3. a) Cycle distribution for 31 non-interacting bosons in a 3D harmonic trap; 
Symbols with error bars denote PIMC results. Analytical results, Eq. (1.36), are 
shown with lines. Cq denotes the permutations of length q. Right panel: Superfluid 
fraction (b) and condensate fraction (c) for varying numbers of 2D trapped bosons. 
Results for boltzmannons ("boltz", B) do not depend on particle number A'^. 

cycles of length 1 occur at higher temperatures (uppermost curve Ci). In this 
regime Bose statistics are of no importance and a PIMC simulation without 
permutation sampling results in the same values for thermal quantities. Note, 
that this does not necessarily imply that the system can be treated classically. 
It rather means that the single-particle wave- functions have a reasonably small 
overlap, while still extending over a finite space. Such quantum particles are 
called "Boltzmannons" . 

The effect of Bose statistics on thermal averages at low temperatures can 
be further investigated by switching it on and off in the PIMC program. In 
the analytical formulas for the noninteracting Bose gas this is easily achieved 
by setting Ci — N and putting every longer cycle number to zero. The impact 
on the superfluid fraction 7s is shown in Fig. 1.3(b). In both cases it reaches 
unity at absolute zero and disappears for high temperatures. The fact that a 
noninteracting system shows a superfluid response may seem surprising since 
superfluidity is a ground state property which appears in systems with a linear 
dispersion and is, thus, inseparably connected with inter-particle interactions. 
However, it does occur even in mesoscopic ideal systems due to the discrete 
nature of the excitation spectrum in finite systems, leading to energy gaps. 

Fig. 1.3(b) indicates that in a finite trapped system superfluidity appears 
also for boltzmannons. Nevertheless, particle exchange allows condensation 
into the ground state at higher temperatures which ultimately increases the 
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Temperature, k^T / h co 

Fig. 1.4. Superfluid fraction for 5 interacting bosons, < A < 100, in a 2D harmonic 
trap versus temperature, T — kBT/tku. Symbols denote PIMC simulations. Dashed- 
doted line displays an analytical result, 7^ = 1 — Iq/Ic, for ideal system, Eqs. (1.40- 
1.41). The inserts show the density distributions at A = 10 and three temperatures: 
denoted as (a),(b),(c). 

superfluid fraction for all temperatures. In contrast to the boltzmannonic 
calculations, the superfluid fraction in the bosonic system, therefore, depends 
on the particle number. An explanation for this behavior can be drawn from 
the analytical expressions for the thermal averages. These depend only on 
single-particle quantities as it should be the case for an ideal system. But an 
exchange cycle of length q contributes to an average like a single particle at a 
g-times lower temperature. The path-integral method uses the same principle 
in the opposite direction — a single-particle is described as a cycle. 

Condensate and superfluid fraction are not identical as it can be seen in 
Fig. 1.3(b),(c). The higher the particle number, the more the condensate frac- 
tion shows its typical polynomial dependence on the temperature (which is 
1 — (T/Tc)^ for 2D trapped systems where Tc is the critical temperature). 
In the thermodynamic limit the "softening" at the transition point vanishes 
completely. The system then has a clearly defined critical temperature be- 
low which Bose-Einstein condensation occurs. On the other hand, taking the 
thermodynamic limit means an infinite system volume which implies hw ^ Q. 
Thus, the superfluid density will disappear completely. The plots in Fig. 1.3 
cannot show this behavior due to the chosen temperature scale which keeps 
hw constant. 

B. Interacting Bosons 
With the PIMC method it is possible to include inter-particle interactions like 
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e.g. Coulomb repulsion (1.24) from first principles. The effective strength of 

the interaction can be controlled by the trap frequency and is measured by 
the parameter A. Fig. 1.4 shows the temperature dependence of the superfluid 
fraction for several values of A (e.g. 2 < A < 10 corresponds to typical particle 
densities in semiconductor quantum dots). The repulsive interaction causes a 
shift of the transition temperatures to lower values. When cooled down, the 
system typically forms a crystal like state in intermediate temperature regions 
until it melts into a ring like structure with delocalized particles, see insets 
in Fig. 1.4. Obviously, the latter shows a high superfluid response which is 
proportional to the ratio of the area enclosed by paths to the cross-section 
of the whole system (sec Eq. 1.26). In the ideal case, the system skips the 
intermediate crystal phase and directly reaches the delocalized state. In the 
case of dominating interaction strengths, the system stays highly localized 
even at absolute zero. Note, that even for the crystal phase the simulations 
yield a non vanishing value 7s . This is a finite size effect because of a nonzero 
area ratio (1.26). 

1.6 Discussion 

a) Let us summarize the main areas of application of the PIMC method: 

• Low temperature systems (relevance of quantum effects) . 

• Small dimensions (system size comparable to De Broglie wavelength A_d). 

• High density: 2-particle distance comparable to or / smaller than Xd- 

b) The simulations allow to solve problem of 1, 2 or A'' quantum particles: 

• Single particle in a complicated potential, e.g. disorder effects effective 
solution of Schrodinger's equation (ground state, T = 0) or finite temper- 
ature extension (density matrix) [7]. 

• Improvement of pair potential at small distances (include quantum ef- 
fects) [9]. 

• Finite number of particles in traps atoms, ions at ultra low temperature 
(Bosc condensates etc.) electrons, holes in quantum dots, e.g. [26]. 

• Macroscopic quantum systems of electrons and ions in astrophysics: planet 
cores, dwarf stars, highly excited semiconductors (many electrons, holes 
in nanostructures), e.g. [15, 16]. 

c) Naturally, this list is far from being complete. It is the generality of this 
first principle approach which makes it increasingly important for many fields 
of physics and beyond. 
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